% This script download magnetic variation data from Super Mag site
% does some quality check, project the data to geomagnetic main field
% vector and calculate the abosulte values


% download with 1) all baselines (daily + annual) removed. 2) only yearly
% removed
% 1) http://supermag.jhuapl.edu/script/download.mag?cli=manojnair+20000101+00%3A00+24%3A00+-s+GDH%2CKUV%2CSKT%2CSTF%2CSVS%2CTHL%2CUMQ%2CUPN+-delta%20none+-baseline+all
% 2) http://supermag.jhuapl.edu/script/download.mag?cli=manojnair+20000101+00%3A00+24%3A00+-s+GDH%2CKUV%2CSKT%2CSTF%2CSVS%2CTHL%2CUMQ%2CUPN+-delta%20none+-baseline+yearly

%2CKOT%2CT41%2CBET%2CT40%2CJCO%2CBET%2CSHU%2CHOM%2CAMU%2CHLM%2CT39%2CTLK%2CKAV%2CARC%2CFYU%2CCMO%2CGAK%2CEAG


%http://supermag.jhuapl.edu/script/download.mag?cli=datat+20031025+00%3A00+240%3A00+-s+ALE%2CKOT%2CT41%2CBET%2CT40%2CJCO%2CBET%2CSHU%2CHOM%2CAMU%2CHLM%2CT39%2CTLK%2CKAV%2CARC%2CFYU%2CCMO%2CGAK%2CEAG+-delta%20none+-baseline+none

%2CDAW%2CINK%2CCPY%2CSAH%2CT08%2CFSP%2CYKC%2CMBC%2CRES%2CCBB%2CALE%2CCY0%2CPGC%2CIQA%2CCDC%2CCHC%2CEKP%2CRAN%2CTAL%2CGHC%2CPBC%2CIGC%2CEKP%2CBLC%2CT35%2CT34%2CCNL%2CT22%2CSMI%2CT37%2CT02%2CFVE%2CT04%2CFMC%2CT36%2CMEA%2CEDM%2CT07%2CRED%2CROT%2CT43%2CT03%2CLET%2CRAL%2CT42%2CSAS%2CT13%2CT05%2CGLN%2CPIN%2CT33%2CFCC%2CT01%2CGIM%2CT11%2CISL%2CT06%2CPIN%2CGLN%2CT32%2CT15%2CT31%2CNAN%2CVLO%2CT30%2CT29
% '+00%3A00+240%3A00+-s-delta%20none+-baseline+none'], ...

%for date_number = datenum(1999,5,10):10:datenum(2012,12,31),
for years = 1995:2012,
    
    date_number = datenum(years,1,1);
    
    date_string = sprintf('%s',datestr(date_number,'yyyymmdd'));
    
    urlwrite(['http://supermag.jhuapl.edu/script/download.mag?cli=manojnair+' date_string ...
        '+00%3A00+240%3A00+-s' ...
        '+ALE%2CDAW%2CINK%2CCPY%2CSAH%2CT08%2CFSP%2CYKC%2CMBC%2CRES%2CCBB%2CALE%2CCY'...
        '0%2CPGC%2CIQA%2CCDC%2CCHC%2CEKP%2CRAN%2CTAL%2CGHC%2CPBC%2CIGC%2CEKP%2CBLC'...
        '%2CT35%2CT34%2CCNL%2CT22%2CSMI%2CT37%2CT02%2CFVE%2CT04%2CFMC%2CT36%2CMEA'...
        '%2CEDM%2CT07%2CRED%2CROT%2CT43%2CT03%2CLET%2CRAL%2CT42%2CSAS%2CT13%2CT05'...
        '%2CGLN%2CPIN%2CT33%2CFCC%2CT01%2CGIM%2CT11%2CISL%2CT06%2CPIN%2CGLN%2CT32'...
        '%2CT15%2CT31%2CNAN%2CVLO%2CT30%2CT292CKOT%2CT41%2CBET%2CT40%2CJCO%2CBET' ...
        '%2CSHU%2CHOM%2CAMU%2CHLM%2CT39%2CTLK%2CKAV%2CARC%2CFYU%2CCMO%2CGAK%2CEAG+'...
        '-delta%20none+-baseline+yearly'], ...
        ['/Users/manojnair/projects/obs_mag_data/Canada/work/original_baseline_removed/' date_string '.txt']);
    fprintf('%s\n',date_string);
    
end;


%% Create station list
% I use awk to select satation
%
% for date_number = datenum(1995,1,1):10:datenum(2012,12,31),

for years = 1995:2012,
    
    date_number = datenum(years,1,1);
    
    date_string = sprintf('%s',datestr(date_number,'yyyymmdd'));
    date_string_year1 = sprintf('%s',datestr(date_number,'yyyy'));
    date_string_year2 = sprintf('%s',datestr(date_number+20,'yyyy'));
    
    system(['awk ''NR <= 62 {next }{if(($1!=' date_string_year1 ' && $1!=' date_string_year2 ') || ($1 == NAN)) {print $1}}''  /Users/manojnair/projects/obs_mag_data/Canada/work/original_baseline_removed/' date_string '.txt' ...
        ' > /Users/manojnair/projects/obs_mag_data/Canada/work/original_baseline_removed/' date_string '_stations.txt']);
    fprintf('Completed station list for %s\n',date_string);
    
end;


% Create an array header  records
% I use awk to select records start with the year of measurements. Print
% the complete time & number of records

% for date_number = datenum(1995,1,01):10:datenum(2012,12,31),
for years = 1995:2012,
    
    date_number = datenum(years,1,1);
    
    date_string = sprintf('%s',datestr(date_number,'yyyymmdd'));
    date_string_year1 = sprintf('%s',datestr(date_number,'yyyy'));
    date_string_year2 = sprintf('%s',datestr(date_number+20,'yyyy'));
    
    system(['awk ''NR <= 60 {next }{if(($1==' date_string_year1 ' || $1==' date_string_year2 ') && ($1 != NAN)){print $1,$2,$3,$4,$5,$6,$7}}''  /Users/manojnair/projects/obs_mag_data/Canada/work/original_baseline_removed/' date_string '.txt' ...
        ' > /Users/manojnair/projects/obs_mag_data/Canada/work/original_baseline_removed/' date_string '_header.txt']);
    fprintf('Completed headers  for %s\n',date_string);
    
end;


% Create data array
% I use awk to select lines starting with anything but an year.
%
% for date_number = datenum(1995,1,1):10:datenum(2012,12,31),

for years = 1995:2012,
    
    date_number = datenum(years,1,1);
    
    date_string = sprintf('%s',datestr(date_number,'yyyymmdd'));
    date_string_year1 = sprintf('%s',datestr(date_number,'yyyy'));
    date_string_year2 = sprintf('%s',datestr(date_number+20,'yyyy'));
    
    system(['awk ''NR <= 62 {next }{if(($1!=' date_string_year1 ' && $1!=' date_string_year2 ') || ($1 == NAN)) {print $2,$3,$4}}''  /Users/manojnair/projects/obs_mag_data/Canada/work/original_baseline_removed/' date_string '.txt' ...
        ' > /Users/manojnair/projects/obs_mag_data/Canada/work/original_baseline_removed/' date_string '_data.txt']);
    fprintf('Completed data file for %s\n',date_string);
    
end;

% if ~all(data(:,7)==data(1,7)),
%     display(date_string);
% end;

%% Now read the data, properly assign the stations to data.


for date_number = datenum(1995,1,1):10:datenum(2012,12,31),
    
    %     for years = 1995:2012,
    %
    %         date_number = datenum(years,1,1);
    
    date_string = sprintf('%s',datestr(date_number,'yyyymmdd'));
    
    % Read the header files
    header_data = load(['/Users/manojnair/projects/obs_mag_data/Canada/work/original_baseline_removed/' date_string '_header.txt']);
    %Read the station list
    
    fid = fopen(['/Users/manojnair/projects/obs_mag_data/Canada/work/original_baseline_removed/' date_string '_stations.txt'],'rt');
    [A] = fscanf(fid,'%s',[3, inf])';
    fclose(fid);
    unique(A,'rows');
    
    
    
    mag_data = load(['/Users/manojnair/projects/obs_mag_data/Canada/work/original_baseline_removed/' date_string '_data.txt']);
    L = mag_data == 999999;
    mag_data(L) = NaN;
    
    
    if (length(mag_data) > 1440 & length(header_data) > 1440),
        
        time_array = datenum(header_data(:,1:6));
        
        if  all(header_data(:,7)==header_data(1,7)) && length(time_array) == length(mag_data)/header_data(1,7),
            
            
            for i = 1:header_data(1,7),
                
                B = A(i:header_data(1,7):end,:);
                
                
                
                data = mag_data(i:header_data(1,7):end,:);
                save(['/Users/manojnair/projects/obs_mag_data/Canada/work/tendaysfiles_baseline_removed/' A(i,:) '_' date_string], 'data', 'time_array');
                
            end;
            
            fprintf('Saved %s\n',date_string);
            
        else,
            
            fprintf('Error %s\n',date_string);
            
        end;
    end;
end;


%% Now make station file for the whole data


S = dir('/Users/manojnair/projects/obs_mag_data/Canada/work/tendaysfiles_baseline_removed/*.mat');
time_array_min = (datenum(1995,1,1,0,0,0): (1/(24*60)): datenum(2012,12,31,23,59,0))';
a = length(time_array_min);




for i = 1:length(S),
    
    filestring(i,:) = S(i).name;
    
end;

[b, ia, ib] = unique(filestring(:,1:3),'rows');%Find the unique stations


for i = 1:size(b,1),
    
    n_data = nan([1,a]);
    e_data = nan([1,a]);
    z_data = nan([1,a]);
    time_array_data = nan([1,a]);
    
    
    L = ib == i; %select a station
    
    filesubset = filestring(L,:);
    
    for j = 1:sum(L),
        
        date_num = datenum(filesubset(j,5:end-4),'yyyymmdd');
        
        copy_index = (date_num - datenum(1995,1,1)) * 1440 + 1;
        
        load(['/Users/manojnair/projects/obs_mag_data/Canada/work/tendaysfiles_baseline_removed/' filesubset(j,:)]);
        
        n_length = length(data);
        
        if n_length > 14400,
            n_length = 14400;
        end;
        
        n_data(copy_index:copy_index + n_length - 1) = data(1:n_length,1);
        e_data(copy_index:copy_index + n_length - 1) = data(1:n_length,2);
        z_data(copy_index:copy_index + n_length - 1) = data(1:n_length,3);
        time_array_data(copy_index:copy_index + n_length - 1) = time_array(1:n_length);
        
    end;
    
    
    
    save(['/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/' ...
        filesubset(j,1:3)], 'n_data', 'e_data', 'z_data', 'time_array_data');
    
    fprintf('Saved %s\n',filesubset(j,1:3));
    
    
    
end;


%% plot ; after converting variation data to absolute by projecting the variation vector to main field (IGRF) vector

[master_a,master_b] = xlsread('/Users/manojnair/projects/obs_mag_data/Canada_array.xlsx');

S = dir('/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/*.mat');
% S = dir('/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files/*.mat');

time_array_min = (datenum(1995,1,1,0,0,0): (1/(24*60)): datenum(2012,12,31,23,59,0))';
a = length(time_array_min);

master_b(1,:) = [];

station_list = cell2mat(master_b(:,2));


for i = 1:length(S),
    
    filestring(i,:) = S(i).name;
    
end;

[b, ia, ib] = intersect(station_list(:,2:4),filestring(:,1:3),'rows');


for i = 1:length(filestring),
    
    load(['/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/' filestring(i,:)]);
    % load(['/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files/' filestring(i,:)]);
    % Calculate IGRF
    icount = 1;
    
    for decyear = 1995:0.1:2013,
        [mag_igrf(icount,:), H, DEC(icount)] = igrf11magm(0,master_a(ia(i),1),master_a(ia(i),2),decyear);
        icount = icount + 1;
    end;
    
    
    % get declination at all points through  interpolation
    declination = interp1(datenum(1995:0.1:2013,1,1), DEC, time_array_min, 'linear');
    cos_dec = cos (declination * (pi/180))';
    sin_dec = sin (declination * (pi/180))';
    
    L = n_data == 0;
    n_data(L) = NaN;
    
    L = z_data == 0;
    z_data(L) = NaN;
    
    
    subplot(311);
    plot(time_array_min,n_data .* cos_dec - e_data .*  sin_dec);
    datetick('x');
    ylabel('X (nT)');
    hold on;
    %     plot(datenum(1995:2013,1,1),mag_igrf(1:10:end,1),'c.');
    title(filestring(i,1:3));
    subplot(312);
    plot(time_array_min,e_data .* cos_dec + n_data .* sin_dec);
    datetick('x');
    ylabel('Y (nT)');
    hold on;
    %     plot(datenum(1995:2013,1,1),mag_igrf(1:10:end,2),'c.');
    subplot(313);
    plot(time_array_min,z_data);
    datetick('x');
    ylabel('Z (nT)');
    hold on;
    %     plot(datenum(1995:2013,1,1),mag_igrf(1:10:end,3),'c.');
    
    
    saveas(gcf,['/Users/manojnair/projects/obs_mag_data/Canada/plots/' filestring(i,1:3) '_Secular_Removed'],'png');
    
    close all;
    
end;

%% Now create netCDF files for these stations
% after converting variation data to absolute by projecting the variation vector to main field (IGRF) vector


[master_a,master_b] = xlsread('/Users/manojnair/projects/obs_mag_data/Canada_array.xlsx');
S = dir('/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/*.mat');
%S = dir('/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files/*.mat');

master_b(1,:) = [];
station_list = cell2mat(master_b(:,1));


for i = 1:length(S),
    
    filestring(i,:) = S(i).name;
    
end;

[b, ia, ib] = intersect(station_list,filestring(:,1:3),'rows');

% get the total stations

localdirectory = '/Users/manojnair/projects/obs_mag_data/Canada/Secular_Removed_netCDF/';

time_array_min = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';
a = length(time_array_min);

for i = 1:length(b),
    
    
    metadatacellarray = master_b(ia(i),:);
    lat = master_a(ia(i) ,1);
    long = master_a(ia(i),2);
    alt = master_a(ia(i) ,3);
    fprintf('Creating %s \n',cell2mat(master_b(ia(i),1)));
    [ncid, X_ID, Y_ID, Z_ID] = Create_Geomag_netCDF(localdirectory, lat, long, alt, metadatacellarray);
    eval(['load /Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/' cell2mat(master_b(ia(i) ,1)) '.mat' ' n_data e_data z_data']);
    % Calculate IGRF
    icount = 1;
    
    for decyear = 1995:0.1:2013,
        [mag_igrf(icount,:), H, DEC(icount)] = igrf11magm(0,master_a(ia(i),1),master_a(ia(i),2),decyear);
        icount = icount + 1;
    end;
    
    
    % get declination at all points through  interpolation
    declination = interp1(datenum(1995:0.1:2013,1,1), DEC, time_array_min, 'linear');
    cos_dec = cos (declination * (pi/180))';
    sin_dec = sin (declination * (pi/180))';
    
    L = n_data == 0;
    n_data(L) = NaN;
    
    L = z_data == 0;
    z_data(L) = NaN;
    
    x_data = n_data .* cos_dec - e_data .*  sin_dec;
    y_data = e_data .* cos_dec + n_data .*  sin_dec;
    
    x_data(isnan(x_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
    y_data(isnan(y_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
    z_data(isnan(z_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
    
    %multiplied by 10 to preserve the decimal points
    
    netcdf.putVar(ncid, X_ID, 0, a,  int32 (x_data .* 10) );
    netcdf.putVar(ncid, Y_ID, 0, a,  int32 (y_data .* 10) );
    netcdf.putVar(ncid, Z_ID, 0, a,  int32 (z_data .* 10) );
    netcdf.close(ncid);
end;

%% plot a map of the stations


[master_a,master_b] = xlsread('/Users/manojnair/projects/obs_mag_data/Canada_array.xlsx');

%S = dir('/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/*.mat');
S = dir('/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files/*.mat');

for i = 1:length(S),
    
    filestring(i,:) = S(i).name;
    
end;

master_b(1,:) = [];

station_list = cell2mat(master_b(:,2));
[b, ia, ib] = intersect(station_list(:,2:4),filestring(:,1:3),'rows');

worldmap([40,85],[-180,-50])
plotm(lat,long)
plotm(master_a(ia,1),master_a(ia,2),'r.','MarkerSize',20);

for i = 1:length(S),
    
    stations(i,:) = ['  ' S(i).name(1:3)];
    
end;

textm(master_a(ia,1),master_a(ia,2), stations);

% Save the station data to a file
fid = fopen('/Users/manojnair/projects/obs_mag_data/Canada/station_list.csv','wt');

fprintf(fid,'Station , Latitude, Longitude \n');
for i = 1: length(b),
    
    fprintf(fid,'%s , %6.2f , %6.2f\n', filestring(i,1:3), master_a(ia(i),1),master_a(ia(i),2));
end;

%% plot for Jasper

nc_fname = '/Users/manojnair/projects/obs_mag_data/YKC_1995_2012.nc';
%nc_fname = '/Users/manojnair/projects/obs_mag_data/YKC_1995_2012_Secular_Removed.nc';
time_array_min = (datenum(1995,1,1,0,0,0): (1/(24*60)): datenum(2012,12,31,23,59,0))';
a = length(time_array_min);

[XYZ, H, DEC] = igrf11magm(0,62+29/60,(245+31/60)-360,2004.5)
[x_data, y_data, z_data , X_ID, Y_ID, Z_ID, ncid] = read_geomag_netcdf(nc_fname, 0, a, 0);
z_data_o = z_data;
load('/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files/YKC.mat');
%load('/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/YKC.mat');
L = time_array_min > datenum(2004,1,1) & time_array_min < datenum(2004,1,2);
ykc_n = n_data(L);
ykc_e = e_data(L);
ykc_t = time_array_min(L);
ykc_x = x_data(L);
ykc_y = y_data(L);
plot(ykc_t, ykc_n);
ykc_z = z_data(L);
ykc_z_o = z_data_o(L);


close all
subplot(311);
plot(ykc_t, ykc_x,'r','LineWidth',2);
set(gca,'FontSize',16);
hold on
plot(ykc_t, ykc_n,'b','LineWidth',2)
datetick('x');
legend('Intermagnet X','SuperMag N');
title(datestr(datenum(2004,1,1)));
ylabel('X/N (nT)');
subplot(312);
plot(ykc_t, ykc_y,'r','LineWidth',2);
set(gca,'FontSize',16);
hold on
plot(ykc_t, ykc_e,'b','LineWidth',2);
legend('Intermagnet Y','SuperMag E');
datetick('x');
ylabel('Y/E (nT)');
subplot(313);
plot(ykc_t, ykc_z_o,'r','LineWidth',2);
set(gca,'FontSize',16);
hold on;
plot(ykc_t, ykc_z,'b','LineWidth',2)
legend('Intermagnet Z','SuperMag Z')
datetick('x');
ylabel('Z (nT)');


% plot rotated data

subplot(211);
plot(ykc_t, ykc_x,'r','LineWidth',2);
hold on;
set(gca,'FontSize',16);
title(datestr(datenum(2004,1,1)));

plot(ykc_t, ykc_n * cos (DEC * (pi/180)) - ykc_e * sin (DEC * (pi/180)),'b','LineWidth',2);
datetick('x');
ylabel('X (nT)');
legend('Intermagnet X','SuperMag X');
subplot(212);
plot(ykc_t, ykc_y,'r','LineWidth',2);
hold on;
set(gca,'FontSize',16);
plot(ykc_t, ykc_e * cos (DEC * (pi/180)) + ykc_n * sin (DEC * (pi/180)),'b','LineWidth',2);
ylabel('Y (nT)');
legend('Intermagnet Y','SuperMag Y')
datetick('x');

% print the data

fid = fopen('/Users/manojnair/projects/obs_mag_data/Canada/YKC.txt','wt');

for i = 1:length(ykc_t),
    
    fprintf(fid,'%s %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f\n', datestr(ykc_t(i)), ykc_n(i),ykc_e(i),ykc_z(i),ykc_x(i),ykc_y(i),ykc_z_o(i));
end;

fclose all;

%% Now test read some netCDF files

time_array_min = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';
a = length(time_array_min);

%ncfname = '/Users/manojnair/projects/obs_mag_data/Canada/Absolute_netCDF/T42_1995_2012_Absolute.nc';
ncfname = '/Users/manojnair/projects/obs_mag_data/Canada/Secular_Removed_netCDF/T42_1995_2012_Secular_Removed.nc';


[x_data, y_data, z_data , dummy, dummy, dummy, dummy] = read_geomag_netcdf(ncfname, 0, a, 0);

subplot(311);
plot(time_array_min, x_data);
axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
datetick('x','keeplimits');
subplot(312);
plot(time_array_min, y_data);
axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
datetick('x','keeplimits');
subplot(313);
plot(time_array_min, z_data);
axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
datetick('x','keeplimits');





